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The  long-term  objective  of  the  research  program  of  which  this  project  was  a  part  is  to  provide  prescrip¬ 
tive  guidance  for  generating  meshes  for  full  aircraft  configurations  that  will  reduce  solution  error  for  a  given 
number  of  degrees  of  freedom  compared  with  meshes  generated  using  today’s  best  practices.  To  achieve  this 
goal,  we  must  be  able  to  identify  features  in  unstructured  meshes  that  have  an  adverse  effect  on  solution 
quality  for  commonly  used  numerical  methods  for  computational  aerodynamics  and  we  must  develop  tech¬ 
niques  to  modify  meshes  to  reduce  those  effects.  During  the  term  of  this  sponsored  project,  we  proposed  to 
take  the  first  step:  identifying  mesh  features  that  negatively  impact  the  accuracy  of  common  discretization 
schemes,  including  identifying  which  discretization  approaches  are  particularly  sensitive  to  variations  in  mesh 
properties. 

This  report  summarizes  the  work  completed  during  the  term  of  this  award.  Full  details  can  be  found  in 
the  publications  cited;  any  that  are  not  readily  available  online  can  be  obtained  by  contacting  the  PI. 


1  Truncation  Error  Analysis 


To  quantify  the  interactions  between  mesh  and  solution  scheme,  we  relied  on  analytic  accuracy  analysis  of  the 
truncation  error,  which  is  defined  as  the  difference  between  the  continuous  partial  differential  operator  and  its 
discrete  equivalent.  Analysis  of  truncation  error  for  structured  meshes  is  based  on  Taylor  series  analysis.  This 
analysis  approach,  which  is  applicable  to  finite  difference,  finite  volume,  and  finite  element  methods  alike,  is 
based  on  the  underlying  assumption  that  the  solution  is  smooth  and  therefore  well  represented  locally  by  a 
Taylor  series  expansion.  This  assumption  is  a  property  of  the  solution,  not  the  mesh,  and  so  the  approach 
should  extend  to  unstructured  meshes.  This  section  gives  a  brief  outline  of  how  this  can  be  done  for  finite 
volume  methods.  For  simplicity  and  compactness,  discussion  here  will  focus  on  linear  model  problems  in  two 
dimensions,  which  can  be  written  as 


du 

dt 


+  C(U)  =  0, 


where  C  ( V )  is  some  linear  spatial  operator  that  an  be  written  in  flux  divergence  form:  C  ( V )  =  V  •  F.  This 
PDE  can  be  discretized  in  space  using  the  finite  volume  method  and  written  as 


dUj 

dt 


Ai 


F  •  ndl 


(i) 


where  Vi  is  the  control  volume  average  of  the  solution  in  control  volume  i,  and  the  right-hand  side  of  the 
equation  is  the  discrete  flux  integral  around  the  two-dimensional  control  volume.  The  computation  of  the 
flux  integral  is  the  key  step  in  determining  spatial  accuracy  for  finite  volume  methods;  regardless  of  the 
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order  of  accuracy  of  the  scheme,  the  flux  integration  process  produces  a  linear  combination  of  control  volume 
averages: 

^  =  -Y/aijUj  (2) 

3 

The  coefficients  aij  are  the  entries  of  the  global  flux  Jacobian  matrix  and  will  depend  on  both  the  problem 
and  the  discretization  scheme.  However  —  for  a  linear  problem  —  the  aij  are  strictly  geometric  terms 
and  therefore  depend  only  the  mesh  topology  and  geometry  in  the  neighborhood  of  control  volume  i.  In  a 
structured  mesh  context,  this  sum  would  be  readily  recognizable  as  a  collection  of  discrete  spatial  derivatives. 
For  unstructured  meshes,  the  sum  is  again  a  collection  of  discrete  spatial  derivatives,  but  much  less  easily 
recognizable,  especially  for  general  stencils  or  for  higher-order  schemes. 

To  determine  the  spatial  accuracy  of  the  scheme,  we  must  evaluate  how  accurately  the  sum  in  Equation  2 
approximates  the  control  volume  average  of  the  linear  operator  C(U)  over  control  volume  i : 

TEj  =  T  ffc  (U)  dA  -  ]T  aijUj .  (3) 


As  written,  the  truncation  error  is  a  combination  of  differential  and  discrete  terms.  These  can  both  be  written 
in  terms  of  derivatives  of  the  solution  at  the  reference  location  Xi  for  control  volume  i.1 
To  begin  this  process,  we  write  the  control  volume  average  ^  of  a  smooth  function  (j)  as 


<t>i  =  T  J  <t>{x  -  Xi)  dA  = 


+  di 


xij  +  TV 


d(j) 

dy 
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Qd 
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where  the  geometric  terms  in  this  equation  are  of  the  general  form 

X™ iTij  =  J-  J  {{x  —  Xj)  +  (xj  —  Xi))n  ■  ((y  —  Uj)  +  (yj  —  Vi))7 
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3 
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(«-*)! 


y  {xj  -  Xi)k  ■  (yj  -  yi)1  ■  xn  kym  lj 


This  general  form  can  be  used  directly  to  replace  the  control  volume  averages  Uj  in  Equation  3  with  a  linear 
combination  of  derivatives  of  the  solution  at  Xi.  For  the  control  volume  average  of  the  linear  operator,  we 
simplify  the  general  case  to  give  a  control  volume  average  in  i: 


d(j) 

dy 


Vi 


dx 2 


x2i  d2(j) 

_  d2(j> 

i  2  +  dxdy 

XVi  +  "o  2 
i  dy2 

(5) 


Equation  5  can  be  applied  to  each  of  the  terms  in  C  (U)  to  evaluate  the  integral  term  in  the  truncation  error. 


1.1  Interior  Schemes 

We  have  applied  this  approach  to  study  the  interaction  between  mesh  quality  and  solution  discretization 
schemes  for  several  model  problems.  Our  initial  work  focused  on  diffusion  operators,  because  of  the  wide 
range  of  schemes  in  common  use  for  such  problems  [6,  5].  We  used  the  mesh  fragments  shown  in  Figure  1 
to  test  the  sensitivity  of  different  discretization  schemes  to  deviations  from  a  perfect,  uniform,  equilateral 
triangular  mesh,  using  cell-centered  control  volumes.  Our  results  showed  that  least-squares  reconstruction 
produces  more  accurate  solution  gradients  —  and  therefore  more  accurate  diffusive  fluxes  —  compared  with 

1  Typically,  the  vertex  for  a  vertex-centered  scheme  or  the  cell  centroid  for  a  cell-centered  scheme. 
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(d)  Stretched  mesh 


(e)  Curved  mesh 


(c)  Scaled  mesh 


(f)  Randomly  perturbed  mesh 


Figure  1:  Uniform  and  distorted  meshes  for  analytic  tests. 


Green-Gauss  gradient  calculations  on  meshes  with  any  nonuniformity.  Also,  we  found  that,  when  using 
second-order  linear  reconstruction,  flux  integrals  were  second  order  accurate  only  on  uniform  meshes.  For 
sheared  and  scaled  meshes,  the  flux  integrals  were  first-order  accurate,  and  for  our  stretched,  curved,  and 
randomly  perturbed  test  meshes,  the  flux  integrals  were  zero  order  accurate.  These  results  are  consistent 
with  other  work  in  the  literature  (for  instance,  [2,  3])  and  also  with  our  subsequent  numerical  experiments 
on  full-scale  unstructured  meshes.  We  also  found  that  averaging  the  cell  gradients  based  on  cell  volume 
when  computing  a  single  face  gradient  for  the  flux  is  consistently  significantly  less  accurate  than  using  either 
arithmetic  averaging  or  weighting  each  cell’s  gradient  with  the  other  cell’s  size.2  Including  Nishikawa’s  jump 
term  [10],  designed  to  account  for  the  difference  in  the  two  reconstructed  solution  values  at  the  face,  in  the 
flux  also  improves  accuracy,  regardless  of  mesh  type.  Unsurprisingly,  we  found  that  the  optimal  value  for  the 
free  coefficient  in  the  jump  term  formulation  is  different  for  unstructured  meshes  than  the  value  Nishikawa 
recommends  for  sturctured  meshes.  This  work  was  later  extended  to  vertex  centered  control  volumes,  with 
similar  results  [17]. 

We  extended  this  analysis  to  advection  problems,  both  the  linear  advection  equation  and  Burgers’  equa¬ 
tion  [5,  7].  Because  Burgers’  equation  is  nonlinear,  the  truncation  error  includes  not  only  terms  like  uXiXj 
but  also  its  cousins  uXiuXj  and  uuXiXj1  where  the  indices  i  and  j  indicate  coordinate  direction.  For  these 
problems,  our  primary  interest  was  in  comparing  the  accuracy  of  upwind  and  centered  flux  formulations. 
Our  results  showed  that,  for  realistic  unstructured  meshes,  there  is  no  statistically  significant  difference  in 
truncation  error  between  these  schemes. 

For  non-linear  problems  more  complex  that  Burgers’  equation,  the  task  of  identifying  all  the  terms  in  the 
truncation  error  expansion  and  sorting  out  their  effects  in  numerical  experiments  becomes  infeasible.  Rather 
than  wasting  effort  by  continuing  to  try  to  solve  an  intractable  problem,  we  turned  our  attention  instead  to 
using  our  truncation  error  estimates  to  improve  solution  error,  as  I  will  describe  in  Section  2. 

1.2  Boundary  Condition  Enforcement 

Our  experience  showed  us  that  the  largest  truncation  error  typically  occurs  at  or  near  the  domain  boundary 
because  of  asymmetry  in  stencils  there,  so  we  turned  our  attention  to  error  analysis  near  boundaries  and  on 

2 The  latter  is  equivalent  to  linear  interpolation  of  the  gradient  in  the  direction  perpendicular  to  the  face. 
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(a)  Solution,  structured  mesh 
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(d)  Solution,  unstructured  mesh 
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(b)  Solution  error,  structured  mesh 


(e)  Solution  error,  unstructured  mesh 


(f)  Truncation  error,  unstructured 
mesh 


Figure  2:  Error  comparison  for  structured  versus  unstructured  meshes  for  a  simple  Poisson  test  case. 


mixed  element  meshes  [12,  11].  We  found  that  the  truncation  error  near  the  boundary  is  strongly  dependent 
on  the  boundary  condition  implementation.  Our  tests  used  a  manufactured  solution,  projected  onto  the 
mesh  by  control  volume  averaging,  and  Dirichlet  boundary  conditions.  The  weakest  boundary  condition 
implementation  we  used  simply  performed  least-squares  reconstruction  as  normal  in  boundary  cells  and 
computed  a  flux  based  on  the  gradient.  This  scheme  uses  no  direct  information  about  the  boundary  condition, 
and  performs  poorly  even  for  this  simple  flux  integral  test.  Adding  constraints  to  enforce  the  boundary 
condition  at  flux  quadrature  points  along  the  boundary  gives  a  strong  boundary  condition  and  reduces 
truncation  error  by  about  a  factor  of  two.  Better  still  is  to  return  to  the  weak  implementation  but  add  to 
the  boundary  flux  a  Nishikawa-style  jump  term  that  accounts  for  the  difference  between  the  reconstructed 
solution  and  the  physical  boundary  condition;  using  the  same  jump  term  coefficient  as  in  the  interior  of 
the  domain  reduces  error  by  about  a  factor  of  four  compared  with  the  strong  BC,  while  a  separate  jump 
coefficient  for  only  the  boundary  faces  can  be  tuned  to  reduce  error  by  a  further  factor  of  about  1.5.  These 
trends  also  apply  at  the  interface  between  quadrilateral  and  triangular  mesh  regions. 

1.3  Eigenanalysis  of  Error 

Error  analysis  for  unstructured  mesh  finite  volume  methods  reveals  two  key  differences  between  error  on 
structured  versus  unstructured  meshes.  First,  for  a  structured  mesh  solution  with  error  of  order  p,  the  trun¬ 
cation  error  is  also  order  p,  whereas  a  unstructured  mesh  solution  of  the  same  order  will  have  truncation 
error  of  order  p  —  d,  where  d  is  the  highest  number  of  derivatives  in  the  governing  partial  differential  equation. 
Second,  the  structured  mesh  truncation  error  is  smooth,  whereas  the  unstructured  mesh  truncation  error  is 
extremely  noisy.  Both  of  these  features  are  illustrated  in  Figure  2.  The  reason  that  unstructured  mesh  trun¬ 
cation  error  is  asymptotically  larger  than  the  discretization  error  is  relatively  well-known:  mesh  irregularities 
prevent  any  cancellation  of  error  from  occurring.  When  this  project  began,  to  our  knowledge,  there  was  no 
rigorous  explanation  of  why  the  discretization  error  is  asymptotically  smaller  than  the  truncation  error.  We 
were  able  to  accomplish  this  by  applying  eigenanalysis  [14,  8,  9]  to  the  generalized  error  transport  equation 
in  discrete  form  [13]: 
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Figure  3:  Truncation  error  for  a  Poisson  problem  in  a  square.  Top  left:  truncation  error.  Top  right:  magnitude 
of  coefficients  in  the  eigendecomposition.  Bottom  row:  eigenvectors  associated  with  the  eigenvalues  labeled 
above. 


where  r  is  the  truncation  error,  5  is  the  discretization  error,  and  is  the  global  flux  Jacobian.  An  example 
eigendecomposition  for  the  truncation  error  for  a  second-order  discretization  for  a  Poisson  problem  is  shown 
in  Figure  3. 

If  we  expand  substitute  eigenvector  decompositions  of  r  =  JT  and  e  =  JT  biXi,  where  the  Xi  are 
right  eigenvectors  of  we  get 
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i 

—  ^  ^ 
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Y  bjXjXj 
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In  other  words,  the  coefficients  in  the  eigendecomposition  of  the  discretization  error  are  smaller  than  their 
counterparts  for  the  truncation  error  by  a  factor  of  the  eigenvalue.  This  result  holds,  of  course,  for  both 
structured  and  unstructured  meshes.  Unstructured  mesh  discretizations  see  an  improvement  in  accuracy 
between  the  truncation  error  and  discretization  error  because  the  dominant  truncation  error  modes  are 
rough  modes,  associated  with  eigenvalues  that  scale  like  l/hd.  For  structured  meshes,  on  the  other  hand, 
the  dominant  terms  in  the  truncation  error  are  smooth  modes,  with  eigenvalues  that  are  0(1).  We  have  also 
been  able  to  use  eigenanalysis  to  explain  why,  in  some  cases,  a  scheme  with  lower  truncation  error  produces 
solutions  with  higher  discretization  error  [8].  As  a  spin-off  from  this  eigenanalysis  work,  we  are  now  beginning 
to  study  mesh  impacts  on  solver  stability  by  examining  the  behavior  of  the  eigenvalues  and  attempting  to 
systematically  modify  the  mesh  to  improve  stability. 
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2  Error  Reduction  for  Unstructured  Mesh  Methods  via  the  Error 
Transport  Equation 

Ultimately,  the  goal  of  CFD  is  to  compute  quantities  of  engineering  interest  accurately  and  efficiently.  In 
computational  aerodynamics,  those  quantities  are  typically  integrals  of  boundary  tractions,  like  lift,  drag, 
yawing  moment,  etc.  The  current  state-of-the-art  for  improving  the  accuracy  of  these  schemes  on  a  single 
mesh  is  adjoint  methods,  which  are  specific  to  the  output  quantity  desired.  There  are  well-known  techniques 
for  combining  the  truncation  error  for  the  primal  problem  with  the  solution  to  the  adjoint  problem  to  improve 
the  accuracy  of  output  quantities.  This  approach  has  been  shown  by  many  researchers  to  be  quite  successful 
for  both  steady  and  unsteady  problems.  For  unsteady  problems,  the  drawback  to  adjoint  methods  is  that 
the  adjoint  problem  must  be  integrated  backwards  in  time,  and  so  the  primal  solution  must  be  available  at 
all  times  to  be  able  to  solve  the  unsteady  adjoint  problem. 

The  twin  disadvantages  of  the  adjoint  method,  then,  are  its  specificity  to  a  single  output  quantity  and 
the  data  storage  challenge  for  the  unsteady  adjoint  problem.  To  address  both  of  these  issues,  we  have  begun 
work  on  solving  the  error  transport  equation,  which  computes  the  discretization  error  in  the  original  physical 
problem.  This  approach  allows  us  to  improve  the  solution  everywhere  (freeing  us  to  compute  any  output 
quantity  more  accurately)  and  is  integrated  forwards  in  time  (making  the  memory  requirement  much  smaller) . 
We  begin  with  a  PDE  written  in  flux  divergence  form: 

dtU  I  V  •/(/’)-  5  (6) 

We  would  like  to  estimate  the  discretization  error  e(x,t)  =  U(x,t)  —  U(x,t)  where  U  is  the  exact  solution  of 
the  PDE  and  U  is  some  approximate  solution.  Substituting  into  the  PDE  and  re-arranging,  we  get  the  error 
transport  equation: 


dt£  +  V  •  (/(e  +  U)~  f(U))  =  -(< %U  +  V  •  fill)  -  S)  (7) 

Note  that  the  source  term  in  this  equation  for  the  discretization  error  is  the  truncation  error  present  for  the 
approximate  solution  U.  If  /  is  nonlinear,  the  difference  on  the  left-hand  side  can  not  be  simplified,  although 
it  can  be  linearized.  In  solving  the  error  transport  equation  discretely,  there  are  three  choices  one  must  make: 
the  order  of  accuracy  p  for  the  primal  flow  problem  (Eq.  6) ;  the  order  of  accuracy  q  for  the  error  transport 
equation  (Eq.  7);  and  the  order  of  accuracy  r  for  the  residual  evaluation  (right-hand  side  of  Eq.  7). 

For  structured  meshes,  there  is  a  relatively  large  amount  of  work  on  the  error  transport  equation;  for 
a  particularly  good  treatment,  see  Banks  et  al  [1].  The  principal  result  is  that,  for  problems  with  smooth 
solutions  and  therefore  well-behaved  truncation  error,  the  order  of  accuracy  of  the  corrected  solution  U  +  s 
(equivalently,  the  order  of  accuracy  of  the  error  estimate  e)  is  min  (p -f  g,  r).  Because  the  computational 
cost  for  the  ETE  is  comparable  to  the  cost  of  the  primal  problem,  this  implies  that,  for  the  cost  of  two 
second-order  flow  solutions  plus  a  single  fourth-order  residual  evaluation,  we  can  get  a  fourth-order  solution. 
In  practice,  for  linear  problems  defect  correction  methods  are  equivalent  to  the  ETE  solution  with  p  =  q. 

The  situation  is  more  difficult  for  unstructured  meshes,  because  the  structured  mesh  behavior  just  dis¬ 
cussed  requires  that  the  solution  be  smooth  and  that  the  truncation  error  be  smooth  and  of  the  same  order 
as  the  discretization  error  [4].  As  we  have  seen,  this  is  not  the  case  for  unstructured  meshes,  so  the  same 
results  can  not  be  expected  to  apply.  Our  experiments  began  with  model  problems  on  randomly  non-uniform 
meshes  in  ID.  We  evaluated  solution  and  gradient  data  for  the  residual  by  using  high-order  reconstruction 
(r  >  p).  To  improve  accuracy  asymptotically,  we  must  found  that  we  had  to  use  q  =  r  >  p,  which  resulted 
in  order  q  accuracy.  In  other  words,  we  had  to  solve  the  error  equation  to  higher  accuracy  than  the  primal 
equation  [16].  This  result  also  holds  for  model  problems  in  multiple  dimensions. 

For  linear  problems,  solving  the  error  transport  equation  to  high  order  is  identical  to  solving  the  primal 
problem  to  high  order,  and  so  is  not  worthwhile.  For  non-linear  problems,  we  have  shown  that  the  steady 
ETE  gains  an  advantage,  in  two  ways.  First,  applying  a  Newton  linearization  of  the  Eq.  7  allows  us,  for 
steady  problems,  to  solve  the  ETE  for  the  cost  of  a  single  linear  solve  of  the  high-order  linearized  system  while 
retaining  the  same  asymptotic  accuracy  as  the  fully  non-linear  solution,  provided  that  2 p  >  g  =  r.  Second, 
the  low-order  primal  solution  is  typically  more  robust  and  faster  to  converge  than  the  high-order  primal 
solution;  the  combination  of  a  low-order  primal  solve  followed  by  the  linearized  ETE  retains  this  robustness 
advantage.  Overall,  we  see  a  speedup  of  a  factor  of  two  or  more  for  using  a  low-order  primal  solution  and 


6 


the  linearized  ETE  to  compute  the  high-order  solution  compared  with  solving  the  high-order  primal  problem 
directly.  These  results  are  based  on  solutions  to  the  Euler  and  laminar  Navier-Stokes  equations  [15] 

We  are  currently  working  to  extend  these  results  to  unsteady  flows.  We  expect  that  we  will  be  able  to 
achieve  structured-like  results  for  the  time  advance,  because  step  size  will  at  worst  vary  smoothly.  We  also 
expect  to  retain  a  robustness  advantage  compared  with  high-order  solutions.  The  comparison  in  effort  and 
results  with  unsteady  adjoint  methods  remains  to  be  seen. 


3  Personnel  Involved 

Several  graduate  students  have  participated  in  this  research.  Alireza  Jalali  (MASc  2012)  did  the  initial  work 
on  error  analysis  for  interior  schemes  on  cell-centered  meshes  (Sec.  1.1)  as  well  as  the  impact  of  eigenstructure 
on  error  (Sec.  1.3). 3  On  the  latter  topic,  he  was  assisted  by  Mahkame  Sharbatdar  (PhD  anticipated  2016). 
Varun  Puneria  (MASc  2015)  performed  the  studies  on  error  near  boundaries  (Sec.  1.2).  Mr.  Gary  Yan,  Mr. 
Puneria,  and  Mr.  Jalali  collaborated  on  the  analysis  of  vertex-centered  methods. 

Gary  Yan  (PhD  anticipated  2017)  began  his  research  work  looking  at  Taylor  series  analysis  for  non-linear 
problems.  Since  that  proved  to  be  a  blind  alley,  he  has  shifted  his  focus  to  the  error  transport  equation. 

Mr.  Puneria  and  Mr.  Yan  have  been  funded  from  this  project  throughout  their  studies  at  UBC.  Ms. 
Sharbatdar  has  been  primarily  funded  from  other  sources,  and  Mr.  Jalali  has  been  entirely  funded  from  other 
sources. 
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